The coupling relationship and driving mechanism between urbanization and ecosystem services in the Yellow River Basin from a multi-spatial scale perspective

Rapid urbanization has led to ecological destruction and associated issues. Macro policies wield substantial influence over urbanization and its relationship with the environment. Without considering the differences in scale, macro policies may be ineffective at addressing urbanization’s adverse impacts on the environment, and even worsen this relationship. We used data on 622 counties, 76 prefectures, and 7 urban agglomerations in the Yellow River Basin to examine the development level, coupling coordination degree, and spatial patterns of urbanization and ecosystem services at three scales during 2000–2020. Further, we explored the driving mechanisms in the relationship between urbanization and ecosystem services. We found that: First, the coupling coordination was relatively low but showed an upward trend. A sizeable spatiotemporal difference existed, with higher (lower) coordination in the east (west). Second, the coupling coordination of each scale exhibited significant spatial positive correlations. The low-value heterogeneous region was embedded around the agglomeration region, and polarization was significant. The larger the scale, the stronger the agglomeration effect. Further, the coupling coordination spatial agglomeration effect of each scale gradually weakened over time. Third, the spatial and temporal distributions of coupling coordination and its agglomeration characteristics at different scales differed. The urban agglomeration scale showed significant overall coordination or agglomeration characteristics, and prefecture and county regions showed local and unique characteristics within urban agglomerations. Fourth, the dominant factors influencing the spatial patterns of the coupling coordination at the county, prefecture, and urban agglomeration scales differed. The interaction and factor detection showed linear and double-factor enhancements. We find that economic development, government policies, environmental protection, and natural factors are the combined effects of urbanization and ecosystem services. Our research method can provide a reference for other river basins, and the results can help governments in formulating policies for sustainable development at different spatial scales.


Introduction
Urban expansion and population agglomeration caused by rapid urbanization cause significant changes in the structure and function of ecosystems [1].Urbanization can bring cities closer, eventually giving rise to urban agglomerations.As the optimal form of spatial organization for urban development, urban agglomerations are a collection of economic agglomerations characterized by rapid industrialization [2].However, the resulting intensification of human activities brought about by the concentration of urban agglomerations has placed tremendous pressures on the ecological environments in urban agglomerations, and those in surrounding counties and cities.These pressures can lead to the loss of arable land [3], reduced carbon stocks [4], and climate change [5], further increasing the risk of ecosystem stress and environmental degradation, and threatening regional sustainable development [6].Essentially, urbanization and ecosystem services exhibit a complex coupling relationship.As such, a key research issue is how should we coordinate this relationship.
In highly urbanized areas, human activities are the dominant factor affecting the functioning of urban ecosystems.Scholars have discussed the relationship between the two at the watershed, province, prefecture, and county scales [7,8], and studied ecosystem services [9], land economic density [10], and human well-being [11] in a single spatial environment.Single-scale research must compare multi-system relationships from multiple spatial perspectives of the same element.
Among studies quantifying the relationship between urbanization and ecosystem services, Fang et al. [12] theorized the interactive coupling effect between urbanization and the ecosystem.Liu et al. [13] proposed the concept of "coupling Rubik's cube" to explain the coupling mechanism between urbanization and ecological environment from the four dimensions of space, time, appearance, and organization.Many studies have quantified the relationship between urbanization and ecosystem services using multi-source data [14], such as yearbook, night light, and ecological index remote sensing data [15].However, there are few regional parameters which quantitatively measure ecosystem service value at the watershed level [16].Here, we use the InVEST model to quantitatively measure ecosystem service capacity.After formulating model coefficients according to regional characteristics, we fit the real data for simulation and use the coupling coordination model to quantitatively evaluate the relationship between the two system elements.
Research on the driving mechanisms of coupling and coordination relationships has mainly adopted methods such as data spatial superposition [17], the STIRPAT model [18], the system dynamics (SD) model [19], the interactive stress model [20], geographical detector [21], and the geographically weighted regression series model [22], among others.These studies show that urbanization and ecological services can interact during its rapid development stage, urbanization significantly impacts the ecological environment.Macro policies to regulate the city can effectively reduce ecological pressures.However, extant research only identifies single factors.Moreover, macro policies do not consider the differences in scale, and consequently, may be ineffective or even worsen the relationship between urbanization and the ecological environment.
In China, the Yellow River Basin (YRB) is an essential experimental area for developing an ecological environment, with ecological protection regarded as an important project.However, balancing urbanization and ecosystem services require better natural conditions and more substantial carrying capacity [23].In 2022, the urbanization rate in the YRB was 61.30%, which is approximately 3.9 percentage points lower than the national average.The proportion of built-up area of core cities in each province is increasing, while the areas of forests and water bodies, and natural ecosystem are gradually shrinking.In particular, in the middle and lower reaches of the YRB, soil erosion and ecological damage are severe, yet urban development continues.As such, the imbalance between urban and ecological development in all provinces and regions along the YRB needs to be solved urgently [24,25].Studies have theoretically discussed the scenario and intensity of the interaction between urbanization and ecosystem services in the YRB to guide regional urbanization and sustainable development of the ecological environment [26].In addition, different urban development policies at different scales have led to deterioration of the ecological environment.Therefore, we need a better explanation of the coupling relationship and change rules between various factors of urbanization and ecosystem services in the YRB from multiple spatial scales.
This study quantifies urbanization and ecosystem services from a multi-spatial scale perspective using data on 622 counties, 76 prefectures, and 7 urban agglomerations in the YRB, and analyzes the coupling coordination degree (CCD) between urbanization and ecosystem services.Further, we study the dominant factors affecting CCD at different spatial scales.Our specific objectives are as follows: (1) To quantify and map the urbanization and ecosystem service levels at the county, prefecture, and urban agglomeration scales in 2000, 2005, 2010, 2015, and 2020 using the coefficient of variation method; (2) Use the coupling coordination model to detect the CCD between urbanization and ecosystem services at the three scales, and use the spatial autocorrelation method explore the spatial agglomeration degree; (3) Identify the dominant factors affecting the CCD at the three scales using factor and interactive detections through the geographic detector model.
This study analyzed and explored the development laws of coupling coordination, explored the influencing factors of coupling coordination from multiple spatial perspectives, and provided decision-making basis for the coordinated development of the two systems.The research helps to provide reference for other watersheds, and the research results also provide a theoretical basis for relevant governments in the Yellow River Basin to formulate urban development and environmental protection policies at different spatial scales.

Overview of the study area
The YRB is an essential ecological reserve and economic zone in China, involving nine provinces (including autonomous regions), including Qinghai, Gansu, Sichuan, Ningxia, Inner Mongolia, Shaanxi, Shanxi, Henan, and Shandong, and seven urban agglomerations, including the Shandong Peninsula, Central Plain, Central Shanxi, Huhhot-Baotou-Ordos-Yulin, Guanzhong Plain, Ningxia, and Lanzhou-Xining urban agglomerations.It includes 76 prefectures and 622 counties (Fig 1 ) and is located at 32˚6 0 53@N~41˚48 0 18@N, 95˚50 0 29@E~119˚06 0 53@E, with a total length of 5464 km and a basin area of 3.6×10 6 km 2 .The average annual precipitation in the whole basin is 466 mm.The area of soil and water loss in the Loess Plateau is 43.4×10 4 km 2 , of which the area with an annual average erosion modulus of more than 5000t/ km 2 is approximately 1.56×10 4 km 2 .The average annual natural runoff of the entire river is 58 billion m 3 , accounting for only 2% of the total river runoff in China [23,25].
Significant economic differences exist in the YRB.In 2021, the GDP of Shandong province exceeded 8 trillion yuan, ranking third in China, while the GDP of Ningxia and Qinghai were only 452.23 and 334.66 billion yuan, respectively [23].The urbanization rate in the YRB (61.30%) is approximately 3.9 percentage points lower than the national average (65.20%),indicating significant regional differences in urbanization.The total land use area of the YRB is 2.03×10 6 km 2 , with grasslands accounting for the most at 41.47% with an area of 8.42×10 5 km 2 .The second common land use is unused land, accounting for 25.86% and an area of 5.25×10 5 km 2 .The YRB comprises complex and diverse ecosystem types, fragile ecological environments, and complex landforms [7].As such, protecting the ecological environment together with rapid urbanization has become essential for high-quality development balanced by ecological protection in the YRB [27].

Data sources
The urbanization rate data of 622 counties, 76 prefectures, and 7 urban agglomerations in the YRB in 2000, 2005, 2010, 2015, and 2020 are extracted from the "Statistical Yearbook," "City Yearbook," and statistical bulletins of provinces and cities in the corresponding years.
Data such as the digital elevation model (resolution of 90 m), land use data (resolution of 30 m), 1 km population density, and GDP density data are from the Resource and Environmental Science Data Center of the Chinese Academy of Sciences.PM2.5 data are collected from the Atmospheric Meteorological Analysis Group of Dalhousie University (http://fizz.phys.dal.ca/).CO 2 data are from the National Geophysical Data Center (http://www.ngdc.noaa.gov/).The administrative boundary map is extracted from the standard map service website of the Ministry of Natural Resources (http://bzdt.ch.mnr.gov.cn/).Road network traffic data are from the public street map platform (http://www.openstreetmap.org/).
Meteorological data, such as precipitation and evapotranspiration, are from the National Meteorological Science Data Center (http://data.cma.cn/).Data on soil type, soil organic matter content, and root depth are from the National Scientific Data Research Center for the Qinghai-Tibet Plateau (http://data.tpdc.ac.cn/).The unified projection coordinates of all raster data are WGS_1984_World_Mercator.According to the InVEST model manual, we export the unified raster data resolution of 300 m to TIF format data and the InVEST model is added to estimate the water production service.This yields the new raster data.Then, the unified zoning statistics are performed to obtain the raster quantification in each county, prefecture, and urban agglomeration unit.

Comprehensive urbanization rate
Following previous studies [11,14], we construct the comprehensive urbanization rate of the YRB from the dimensions of population, economy, land, and society (Table 1).We divide the We use the coefficient of variation method to calculate the index weight [28]: W i is the weight of the i th index, and V i S i � A i W i are the coefficient of variation, standard deviation, and average of the i th index, respectively.

Quantifying ecosystem services
We use the InVEST model to assess ecosystem services in the YRB considering region-specific parameters of the actual ecosystem.Following relevant studies on ecosystem services [29,30], we use five modules to evaluate the ecosystem service capacity of the YRB: carbon sequestration service (Y 1 ), habitat quality (Y 2 ), soil conservation (Y 3 ), water conservation (Y 4 ), and food supply (Y 5 ) (Table 2).Then, following existing studies [31], we use SPSS 27 and the K-means clustering method to classify the ecosystem service capability of the YRB into four categories: primary (0.00, 0.15), intermediate [0.15, 0.39), intermediate, advanced [0.39, 0.53) and advanced [0.53, 1.00).

Carbon sequestration.
The carbon sequestration parameters of different land use types are quantified using the InVEST model carbon module [32]: where C x is the carbon sequestration of grid unit x (t/hm 2 ), C above is the carbon stored in the aboveground biomass, C below is the carbon stored in the underground biomass, C soil is the carbon stored in the soil, and C dead is the carbon stored in dead organic matter.

Habitat quality.
The habitat quality module of the InVEST model is used to quantify the threat and sensitivity characteristics, degree of habitat degradation, and habitat quality under spatial and temporal changes [33].
where H j is the habitat suitability of land use type j, D 2 xj is the total threat level of grid pixel x in land type j, and k is a saturation constant.

Soil conservation.
The soil conservation module of the InVEST model is based on the soil loss equation.It is calculated based on data such as landform, climate, vegetation, and land management status [34] as follows: RKLS is potential soil erosion, RULSE is substantial soil erosion, SC is actual soil conservation (tkm -2 a -1 ), R is the rainfall erosivity factor, and K is the soil erodibility factor.LS is the slope factor, L is the length, and S is the slope.C is the vegetation management factor, and p is the factor of soil and water conservation measures.

Water conservation.
The water yield module of InVEST model quantifies the water yield.According to Budyko theory, the ratio between the actual evaporation and precipitation is related to potential evapotranspiration [35].The annual water yield of each pixel WY (x) on the ground is calculated as follows: where AET(x) is the actual evapotranspiration of x pixel and P(x) is the annual precipitation on x pixel.

Grain supply.
Research shows a significant linear relationship between grain production and the normalized vegetation index.We simulate the service space of grain production in the YRB using statistical data of grain production in the YRB and NDVI data [36] as follows: where CP x is the grain yield of the grid unit (t), NDVI x is the NDVI of the grid unit x, NDVI sum is the normalized difference vegetation index of the entire study area, and CP sum is the total grain yield of the study area.

Coupling coordination degree model
We compute CCD [15,16] as follows: where u 1 and u 2 are the urbanization rate and ecosystem service level, respectively; D is the CCD; and T is the degree of coordination of urbanization and ecosystem services.α and β are the contribution of urbanization and ecosystem services, respectively; we set α = β = 0.5.Following existing research [14,37], we divide the CCD into severe disorder, mild disorder, moderate, coordination, and superior coordination at intervals of 0.4, 0.6, and 0.7 coordination levels.

Bivariate spatial autocorrelation
We used Moran 's I index to measure global spatial autocorrelation [38]: where i is the global autocorrelation index; x i and x j represent the observed values of units i and j, respectively; n is the number of spatial samples; and W ij is the spatial weight matrix.Moran's I = 0 indicates that the spatial distribution pattern tends to be random.Moran's I > 0 indicates that space is positively correlated.The larger the value, the more substantial the spatial similarity.LISA is used to measure the local spatial autocorrelation and divided into four types: high-high (H-H), high-low (H-L), low-high (L-H), and low-low (L-L).

Geodetector model
Factor and interaction detections in the geographic detector model [36] are used to identify the influencing factors and their interactions as follows: where q is the explanatory degree of each driving factor of urbanization; nσ 2 represent the total sample size and variance, respectively; and n h σ h represent the sample size and sample variance of the h-th, respectively.The closer the value of q is to 1, the greater the influence of this factor on multi-scale CCD, and vice versa.
Interaction detection [21] is the explanatory force value obtained by the coupling coordination degree of two systems under multiple spatial perspectives, which is jointly influenced by two different driving factors.The explanatory power of interactive detection is enhanced or weakened, or the two are independent, compared to the explanatory power of factor detection.There are five main categories of interaction factor detection relationships (Table 3).
Following relevant scholars [14,39,40], this study considered the CCD between multiscale urbanization and ecosystem services in the YRB as the explained variable.To establish the index, we select explanatory variables from six perspectives: economic development, government behavior, social investment, essential public services, environmental protection, and natural conditions (Table 4).We then use a geographical detector to detect the factors and Table 3. Types of interactive detection results for geographic detectors.

Theory of Judgment
Interaction interactions of the CCD between the county, prefecture, and urban agglomeration scales in the YRB.

Spatial distribution characteristics of urbanization and ecosystem services in the YRB
The urbanization development level in the YRB is relatively low (   (63.67%).The prefecture scale accounts for 34.21% of the total number of cities, mainly distributed in the Shandong Peninsula and Central Plains urban agglomerations, contrary to the spatial distribution of urbanization levels.This is because the expansion of urban space enhances the intensity of human activities, leading to the destruction of the ecological environment and ecological problems, and reducing the level of ecosystem services.
In the high/low level of overall urbanization and ecosystem service in urban agglomerations, the levels of urbanization and ecosystem services at the prefecture and county scales differ.This indicates that the overall development level should be evaluated from different spatial perspectives.The larger the scale, the more macro the perspective.It does not indicate that the prefecture or county levels are also high CCD at a microscale.This is because the spatial and temporal heterogeneity of the distribution of economics, land, population, and other resources within the prefecture and county, as well as the degree of ecological damage, will lead to changes in the factors used to evaluate urbanization and ecosystem services, thus affecting the overall development level.For example, if the overall urbanization level of the Shandong Peninsula urban agglomeration is relatively high, cities with different levels of urbanization development will exist within it, and the micro-scale in county will once again show different urbanization development levels from the prefecture.

Spatiotemporal pattern of CCD between urbanization and ecosystem services in the YRB
From the perspective of developmental trends (Fig 3A ), the spatiotemporal differences in the CCD at different scales are large and fluctuation characteristics are significant.From 2000 to 2020, the CCD of urbanization and ecosystem services at different scales in the YRB is county (0.37) < prefecture (0.38) < urban agglomeration scales (0.54).The CCD at the county scale increases by 0.07 from 2000 (0.37) to 2020 (0.44) albeit in a U-shaped trend of first decreasing and then increasing.The prefecture scale decreases by 0.03 from 2000 (0.40) to 2020 (0.37).The urban agglomeration scale increases by 0.11 from 2000 (0.50) to 2020 (0.61), showing a "V"-shaped fluctuation.Overall, the CCD between urbanization and ecosystem services at the county and urban agglomeration scales increases in the YRB, peaking in 2020, with an overall stable albeit slight decrease at the prefecture scale.From the location of the numerical concentration of the CCD (box position of box plot), we get: county scale < prefecture scale < urban agglomeration.This indicates that within the same spatial scope, different spatial scales will show different agglomeration characteristics.
From 2000 to 2020, multiple spatial scales were found in the quantity proportion of each coupling coordination stage (Fig 3B).For all the three spatial scales, the number of samples in the moderate coordination stage was the largest, and the quantity proportion showed a trend of decreasing.The number of samples in the moderate coordination stage began to increase in 2010, and the number of samples in the excellent coordination stage accounted for the smallest, which remained basically unchanged.Strong stability.At the county scale, the number of counties in the moderate imbalance stage reached 61.60% (2000), and the number of counties in the moderate imbalance stage decreased to 40.00% in 2020, and the number of counties in the moderate coordination stage increased by 22.72 percentage points from 2000 (10.61%) to 2020 (33.33%).The number of cities in each coupling coordination stage is relatively stable, and the ratio changes little.Due to the relatively large scale and small quantity of urban agglomerations, the number of seriously dysfunctional urban agglomerations has decreased to 0 since 2010, and the proportion of moderately coordinated urban agglomerations has increased by 14.34 percentage points from 2010 (14.33%) to 2020 (28.67%).This indicates that the coupling coordination of urbanization and ecosystem services at different spatial scales in the Yellow River Basin is on the rise from 2000 to 2020, and the smaller the spatial scale, the smaller the mean coupling coordination degree.
Spatially, the distribution of CCD at different scales significantly differs between the east and west, and north and south.Areas with severe county imbalance are concentrated in central Shanxi, Guanzhong Plain, Ningxia along the Yellow River, and the Lanzhou-Xining urban agglomeration, and generally spread in the upper and middle reaches of the Yellow River (  Next, the urbanization and ecosystem service levels of some counties and prefectures scattered in the coastal area of Shandong province are high, exhibiting high level and good coordination.Meanwhile, the areas with severe disorder at the prefecture level show the same spatial distribution characteristics as those at the county level.Moderately coordinated prefectures are concentrated in the Shandong Peninsula, Central Plain, and Huhhot-Baotou-Ordos-Yulin urban agglomerations.Seriously uncoordinated prefectures are concentrated in Ningxia along the Yellow River and Lanzhou Xining urban agglomeration, while moderately uncoordinated prefectures are scattered in the Guanzhong Plain urban agglomeration and concentrated in central Shanxi urban agglomeration.In 2019, the high-quality development plan for the YRB sought to rectify areas with low levels of coordination, paying more attention to ecological protection in the process of urbanization, and abandon polluting and consumptive industries to drive urban economic development.In 2020, the Shandong Peninsula and Central Plains urban agglomerations transformed from moderate to excellent coordination.These areas have always been in excellent coordination with solid stability.This shows that analyzing spatial changes from a macro perspective can better highlight the overall characteristics and changing trends.The spatial distribution characteristics of the CCD at different scales differ.The CCD is more stable at the urban agglomeration scale than that at the prefecture and county scales, and the change is the smallest.However, the smaller the size, the stronger the regional heterogeneity.Notably, the urban agglomeration scale is in the stage of excellent coordination.Meanwhile, the prefecture and county scales exhibit moderate coordination or disorder.This shows that the overall coordination of urbanization and ecosystem services at the macro-scale does not prove that the internal prefecture-county system is also coordinated.Mesoscale prefectures and microscale counties show unique characteristics of urbanization development and ecosystem services.This may be because the macroscale perspective often ignores the existing defects in the urbanization and ecosystem service development process from the mesoscale or microscale perspective.

Spatial correlation characteristics of CCD between urbanization and ecosystem services in the YRB
Next, we measure the global spatial autocorrelation of the CCD between urbanization and ecosystem services from 2000 to 2020 and explore the spatiotemporal heterogeneity of the CCD.The spatial dependence of the CCD at all spatial scales in the YRB is significant, showing a positive spatial correlation.The Moran's I values at multiple scales in 2020 are county (0.413) < prefecture (0.534) < urban agglomeration scales (0.670).The spatial aggregation of the CCD at all scales has weakened over time, indicating that the agglomeration effect of core cities has weakened.
Next, we conduct a LISA agglomeration analysis, dividing spatial agglomeration categories into four regions: H-H, H-L, L-H, and L-L (Fig 5).From 2000 to 2020, all three scales show significant agglomeration changes, obvious polarization, and good spatial stability.The countylevel H-H agglomeration has expanded from the central and eastern regions of the Shandong Peninsula to its whole prefectural group.In 2020, it spread to the urban agglomeration of the Shandong Peninsula, and the central and western regions of the Central Plain urban agglomeration.The H-H agglomeration-type area of Huhhot-Baotou-Ordos-Yulin urban agglomerations has expanded from Hangjin Banner to the county of Ordos.In 2015, Kundulun District and Tumed Right Banner of Baotou City were H-H agglomeration areas.
The L-L agglomeration type is dominant, accounting for 15.59-39.54% of the total number of county units, mainly distributed in the central Shanxi urban agglomeration, central and northern Ningxia urban agglomeration along the Yellow City agglomeration, Guanzhong Plain urban agglomeration, and the eastern counties of Lanzhou-Xining urban agglomeration.It is scattered in the central Shanxi urban agglomeration, northern edge of Guanzhong Plain urban agglomeration, and exhibits dot distribution in the middle of Lanzhou Xining urban agglomeration.
The H-L agglomeration areas are mainly distributed at the edge of the L-L agglomeration area, with slight spatial variation, with their numbers gradually increasing to 87 counties in 2020.From 2000 to 2015, L-L agglomeration prefectures were mainly distributed in areas with good development potential in the middle and lower reaches of the Yellow River.The highquality coupling and coordination areas gradually became consistent with county and urban agglomeration scales.The spatial agglomeration area at the urban agglomeration scale shows little change.The H-H agglomeration area is mainly distributed in the Shandong Peninsula and Central Plains urban agglomerations.The L-L agglomeration area is distributed in central Shanxi and Ningxia along the Yellow River prefecture group, with little spatial change over time.The stability of spatial agglomeration is the strongest at multiple scales.
In summary, by 2020, H-H agglomeration areas at different scales are generally distributed in the Shandong Peninsula and Central Plains urban agglomeration, whereas L-L agglomeration areas are mainly distributed in central Shanxi, the Guanzhong Plain, and local areas of the Lanzhou Xining urban agglomeration.The larger the scale of the study area, the more macro-agglomerations are observed.The non-significant areas at the urban agglomeration scale show significant spatial agglomeration characteristics at the county and prefecture scales.This is an essential geographical basis for future research on targeted and personalized policy suggestions for urban development.

Driving mechanism of spatial pattern of CCD between urbanization and ecosystem services 4.4.1 Selection of influencing factors and identification of leading factors.
The explanatory variables for county, prefecture, and urban agglomeration scales are imported into the geographic detector model.The R program is used to select the optimal discrete classification method, and we obtain the influence q value of the explanatory variables at each scale on the CCD of the two systems in the YRB (Table 5).
The mean explanatory power of the influencing factors across scales is as follows: county (0.11-0.26) < prefecture (0.13-0.55) < urban agglomeration scales (0.09-0.94).Economic development, social investment, and environmental protection factors play a leading role in the CCD of spatial patterns of multiscale urbanization and ecosystem services.Essential public services and government behavior play an important role, while natural factors have the smallest influence.The average slope of the western urban agglomeration in the YRB exceeds 10˚, level of urbanization development is low, ecological damage is minor, and CCD is small.The eastern plain region has concentrated population, promising economic development, a good foundation for urban development and construction, high technology and environmental protection ability, and high coupling and coordination.
The dominant factors at each scale are different.Those at the county scale are GDP per capita, CO 2 emissions, and total social fixed asset investment.The dominant factors at the municipal scale are CO 2 emissions, road length, and the savings balances of urban and rural residents.The dominant factors of the urban agglomeration scale are the balance of urban and rural residents' savings, investment in society's fixed assets, and the proportion of the added value of secondary and tertiary industries in the GDP.
The factor detection results show that, as the spatial scale increases, the explanatory power of GDP and fiscal expenditure per capita gradually increases.In particular, the explanatory power of GDP per capita is ranked first at the county scale, sixth at the prefecture scale, and last at the urban agglomeration scale, indicating that the micro-influencing factors of per capita GDP have greater explanatory power at the micro-scale.The proportion of the added value of secondary and tertiary industries in the GDP, balance of savings of urban and rural residents, and investment in fixed assets of the whole society gradually rank at the top in terms of explanatory power as the scale increases.The industrial structure is an essential macro dimension for measuring economic development.Meanwhile, the balance of savings is an important indicator of the disposable income of residents and people's living standards.Investment in the fixed assets of society further adjusts the economic structure and regional distribution of productive forces by establishing new sectors, enhancing economic strength, and creating material conditions for improving people's material and cultural life.The three economic indicators and macro impact factors of social investment have excellent explanatory power at the urban agglomeration scale.The influencing factors with excellent explanatory power at the municipal scale include road length (national and provincial roads) and CO 2 emissions.The spatial scale of the explanatory power of road length is as follows: prefecture (2 nd ) > urban agglomeration (4 th ) > county scales (7 th ).Dependence on transportation accessibility is significantly lower at the county scale than at the prefecture scale.As an essential factor for measuring environmental protection, CO 2 emissions is a mesoscale factor and its explanatory power at the municipal scale is strong.

Interactive detection results.
Each influencing factor has an interactive effect on the evolution process of the spatial pattern of the CCD between urbanization and ecosystem services at the county, prefecture, and urban agglomeration scales from 2000 to 2020.The explanatory power significantly increases after the pairwise interaction of different influencing factors.The multiscale interaction detection results show different degrees of two-factor or nonlinear enhancement without weakening or independent relationships (Fig 6).
The spatial pattern of the coupling relationship between multiscale urbanization and ecosystem services in the YRB is due to the interactions between economic development, government behavior, social investment, essential public services, environmental protection, and natural factors.
The interactions between CO 2 emissions, GDP per capita, and other factors are strong at the county scale.The explanatory power is more than 44%, indicating that the development of the CCD spatial pattern of urbanization and ecosystem services in the YRB is mainly affected by economic development, environmental protection, and other factors.Next, the influence of the number of primary and secondary school students at the county scale is weak in single-factor detection.In contrast, interactive detection shows a nonlinear enhancement after interaction with other factors, especially when the q values of the interaction with GDP per capita and road length reach 0.849 and 0.682, respectively.Only the education level has a long and slow effect on the spatial pattern of CCD at the county scale.At the prefecture scale, the savings balance of urban and rural residents, road length, and CO 2 emissions strongly interact with other factors, and the explanatory power q value exceeds 70%.The interactive detection of all factors at the prefecture scale shows nonlinear enhancement.The explanatory power of environmental protection, traffic accessibility, and government per capita fiscal expenditure is much higher than that of factor detection, peaking at 0.883.This shows that municipal government behavior, environmental factors, and fundamental public service policies significantly impact other factors, and their leading role is significant.The urban agglomeration scale's q-value in the factor detection stage has the most substantial multiscale explanatory power, and interactive detection has the most interaction types enhanced by double factors.Among them, the proportion of per capita GDP and the added value of secondary and tertiary industries in the total added value, per capita GDP and average slope, and society's fixed asset investment and average slope show nonlinear enhancement.The q values after the interaction are 0.248, 0.510, and 0.889, respectively.This shows that per-capita GDP is at the bottom of the scale factor detection of urban agglomeration, which needs to be fully reflected by the combined effects of macroeconomic influencing factors, social investment, and natural factors.Meanwhile, it reflects that adjusting the social investment level has a certain effect.The transformation of economic structure plays a significant role in improving people's living standards.Economic and ecological sustainable development policies, such as strengthening environmental protection and optimizing the economic structure.This is conducive for improving the CCD between urbanization and ecosystem services in the YRB to excellent coordination.The explanatory power of GDP per capita and government fiscal expenditure per capita at the county scale is substantial, reaching 0.43, and interaction detection with other factors shows nonlinear enhancement.The CCD between urbanization and ecosystem services at the county scale is directly driven by GDP per capita and government behavior.The explanatory power of road length and CO 2 emissions at the municipal scale is strong, reaching the highest value of 0.74.Traffic accessibility and environmental protection are necessary to promote the flow of production factors and enhance the sustainable development of the economy, society, and ecology, directly affecting the CCD of urbanization and ecosystem services at the municipal level.At the urban agglomeration scale, society's fixed asset investments, and the savings balance of urban and rural residents have strong explanatory power.
In conclusion, the dominant factors influencing the spatial pattern of the CCD between multiscale urbanization and ecosystem services in the YRB are different.The dominant influencing factors at the county scale are micro-influencing factors directly related to people's living standards and quality.At the prefectural scale, the meso-level enhances the livability characteristics of the prefecture, which can ensure the overall urban and rural living standards and quality of influencing factors.Thus, considering the economic and industrial structures, and social asset investment, the CCD between urbanization and ecosystem services in the YRB is comprehensively affected by various factors, such as the economy, government policy, and environmental protection, at the micro, meso, and macro levels.

Discussion
Significant spatial and temporal differences exists in the level of multiscale urbanization and ecosystem service capacity in the YRB from 2000 to 2020.This study discusses the spatial distribution characteristics of urbanization and ecosystem services in the YRB at different spatial scales.The city's mesoscale and county's microscale show different system-level distributions from the macroscale of the urban agglomeration.This indicates that the larger the scale, the more macro the research perspective, and the results show strong integrity characteristics.Among the characteristics of the system's overall high/low level, there are heterogeneous characteristics within the prefecture and county.By comparing the development levels of urbanization and ecosystem services from different spatial perspectives, we make reasonable suggestions to compensate for the unbalanced development of county urbanization and protect the ecological environment [41].
Analyzing the spatiotemporal characteristics at different scales, the subsystem development level of the two systems is closely related to the CCD level of the whole system.The larger the scale, the higher the degree of CCD.We find an upward trend in the coupling coordination at various scales.However, the overall level is not high, and mainly manifests as a transition from moderate imbalance to moderate coordination.The intra-regional differences across scales are ordered as follows: urban agglomeration < prefecture < county.This is consistent with prior research [37].We highlight the stages and heterogeneity of the coupling coordination of the two systems in the multi-scale coupling and coordination development of urban agglomerations, prefectures, and counties.Further, we provide theoretical explanations for the severe polarization phenomenon and intensification of regional imbalance.We propose a geographical basis for future research and development of government policies, and fundamental changes in the coupling relationship between urbanization and ecosystem services.
Based on multiscale CCD driving mechanism analysis, Wang and Xu [39] mainly referred to the principle and prospect of geographical detectors.Importantly, we find that the larger the scale, the greater the explanatory power of the influencing factors.However, differences in spatial scale lead to different dominant influencing factors.Per capita economic and financial indicators have a more significant impact on county-scale CCD, while environmental protection and transportation factors have a more significant impact on prefecture-scale coupling and coordination [42].Economic structure, social investment, and other factors have a more significant impact on the urban agglomeration scale.Different economic indicators also have different explanatory power at different spatial scales.This differs from the research results driven by a single spatial scale economic factor [43,44].Further, we find that regions with high economic development levels, reasonable industrial structures, developed infrastructure, and better environmental governance have more resource allocation advantages and promote system coupling and coordination development, consistent with Li and Zhang [16].
This study has some limitations.First, the statistical caliber of the relevant data in the different regions was somewhat consistent.Second, we have only conducted our analysis using an existing index system, which may lead to errors in the conclusions.Third, we have considered each unit of the administrative region within the county, prefecture, and urban agglomeration scales as an independent research object while exploring the level of CCD.Future works can integrate multi-source remote sensing data to enrich the rationality and reliability of the demonstration data.Further, owing to data availability concerns, there is no long time series for continuous years of research; thus, researchers should confuct precise and in-depth research.
At the urban agglomeration scale, national macro policy, fiscal revenue, and expenditure policy, such as gross national product, industrial added value, and other macro indicators, can be used to guide the overall direction of high-quality urban development and ecological protection in the YRB.
At the county level, stakeholders and policymakers can focus on specific micro dimensions, such as per capita GDP, and per capita fiscal revenue and expenditure.For good urban agglomeration development, local underdeveloped counties should undertake targeted development efforts focused on the key influencing factors at the county scale.Further, specific methods of county urbanization and ecological protection can be formulated according to the county's regional characteristics.
At the prefecture scale, spatial policy planning should be conducted from the meso-dimensions of the degree of urban greening, urban construction, urban environmental protection, and overall urban and rural development to create a good link and bridge between the macro perspective of urban agglomerations and the micro perspective of counties.Efforts can improve urban priorities, and promote the high-quality development of urbanization and ecological protection of urban agglomerations in flow regions.

Conclusion
Promoting the coupled and coordinated development of urbanization and ecosystem services is essential for sustainable development.Considering the multi-spatial scales in the YRB, this study uses spatial autocorrelation analysis to explore the spatial evolution characteristics of the coupled and coordinated development of multiscale urbanization and ecosystem services in the YRB from 2000 to 2020.We use geographic detectors to identify the main influencing factors and their interactions.Our findings are outlined below.
First, the spatiotemporal differences in the CCD between urbanization and ecosystem services at different scales are significant and show an upward trend.The larger the scale, the greater the CCD.The overall distribution characteristics of the CCD at the county, prefecture, and urban agglomeration scales are consistent, showing that the east is high and the west is low.However, the county and prefecture regions show noticeable regional particularities: the area with high coupling of the entire urban agglomeration, and the local area with low CCD of the inner prefecture and county of the urban agglomeration.
Second, we find a significant positive spatial correlation between the CCD at each spatial scale.The polarization characteristics of the CCD at different scales are significant with good spatial stability.The larger the scale, the stronger the agglomeration effect.Over time, the spatial agglomeration of the CCD gradually weakens and the agglomeration effect of core cities also weakens.The local characteristics of county-scale agglomerations are opposite to those of prefectures and urban agglomerations.The smaller the scale, the more specific the particularities.
Finally, the influencing factors differ across scales.The larger the scale, the stronger the explanatory power of the influencing factors.The interactive probing results show varying degrees of two-factor or nonlinear enhancement.Natural and socioeconomic factors drive the spatial pattern of the CCD; the natural conditions have gradually weakened and socioeconomic conditions have gradually strengthened.

Fig 1 .
Fig 1.Seven urban agglomerations in the YRB.https://doi.org/10.1371/journal.pone.0293319.g001 Fig 2).The number of counties in the rapid and high-speed development stages is small, and they are scattered in the Shandong Peninsula and Central Plains urban agglomerations.Spatially, the urban and urban agglomeration scales are in the budding stage, with a large number of urban agglomerations concentrated in the Lanzhou-Xining, Guanzhong Plain and central Shanxi, Ningxia along the Yellow River, and Huhhot-Baotou-Ordos-Yulin urban agglomerations in the initial stage, and Central Plains and Shandong Peninsula urban agglomerations in the rapid development stage.The level of ecosystem services is also low.The spatial distribution of ecosystem services is opposite to that of the urbanization development level.The majority of counties are low-level

Fig 2 .
Fig 2. Evaluation of multi-scale urbanization and ecosystem services in the YRB.U-A refers to urban agglomeration.https://doi.org/10.1371/journal.pone.0293319.g002 Fig 4).Most counties in the central Shanxi Province are resource-based and have poor ecological environments.The Ningxia urban agglomeration along the Yellow River and Guanzhong Plain has a poor ecological environment.Due to its high ecosystem service index and low urbanization level, the Lanzhou Xining urban agglomeration suffers from a severe imbalance of the urbanization lag type.The mild disorder areas are mainly distributed at the edges of the severe disorder areas.They are concentrated in the Shandong Peninsula urban agglomeration, Central Plains urban agglomeration, and Huhhot-Baotou-Ordos-Yulin urban agglomerations.Moderately coordinated regions are distributed in the Shandong Peninsula, Central Plain, and Huhhot-Baotou-Ordos-Yulin urban agglomerations, with their number increasing in 2020.Notably, the Qinghai-Tibetan Autonomous Prefecture exhibits excellent

Fig 3 .
Fig 3. CCD proportion diagram of multi-scale urbanization and ecosystem services in the YRB.U-A refers to urban agglomeration.a) Box plot of the temporal distribution of multi-scale CCD.b) Proportion of the number of samples in each stage of multi-scale CCD.https://doi.org/10.1371/journal.pone.0293319.g003

Fig 5 .
Fig 5. Distribution of spatial correlation types of CCD between urbanization and ecosystem services across scales in the YRB.U-A refers to urban agglomeration.https://doi.org/10.1371/journal.pone.0293319.g005

Fig 6 .
Fig 6.Interactive detection results of the driving mechanisms of the CCD between the multi-scale two systems in the YRB from 2000 to 2020.Note: The interaction in Fig 6(A) and 6(B) is nonlinear enhancement.In Fig 6(C), *denotes nonlinear enhancement and the rest is two-factor enhancement.https://doi.org/10.1371/journal.pone.0293319.g006

4 . 4 . 3
Formation mechanism of the coupling coordination relationship between urbanization and ecosystem services.Based on the results of dominant factor identification and interaction detection, we construct the formation mechanism of the spatiotemporal evolution of the CCD between multiscale urbanization and ecosystem services in the YRB (Fig 7).

Table 5 . Results of geographic detector model on the influencing factors of CCD of the two systems at multiple scales in the YRB. Driving factors County q value Prefecture q value Urban agglomeration q value 2000 2005 2010 2015 2020 Order 2000 2005 2010 2015 2020 Order 2000 2005 2010 2015 2020 Order
Note: The Z 11 index at the county scale only passes the test in 2020.The P value at the prefecture scale in 2000 is 0.421, which fails the test. https://doi.org/10.1371/journal.pone.0293319.t005